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Abstract 

The objective of this paper is to define empirical parameters (or closure models) for an 
initial thermal non-equilibrium porous-media model for use in Computational Fluid 
Dynamics (CFD) codes for simulation of Stirling regenerators. The two CFD codes 
currently being used at Glenn Research Center (GRC) for Stirling engine modeling are 
Fluent and CFD-ACE. The porous-media models available in each of these codes are 
equilibrium models, which assume that the solid matrix and the fluid are in thermal 
equilibrium at each spatial location within the porous medium. This is believed to be a poor 
assumption for the oscillating-flow environment within Stirling regenerators; Stirling 1-D 
regenerator models, used in Stirling design, use non-equilibrium regenerator models and 
suggest regenerator matrix and gas average temperatures can differ by several degrees at a 
given axial location and time during the cycle. A NASA regenerator research grant has been 
providing experimental and computational results to support definition of various empirical 
coefficients needed in defining a non-equilibrium, macroscopic, porous-media model (i.e., to 
define “closure” relations). The grant effort is being led by Cleveland State University, with 
subcontractor assistance from the University of Minnesota, Gedeon Associates, and 
Sunpower, Inc. Friction-factor and heat-transfer correlations based on data taken with the 
NASA/Sunpower oscillating-flow test rig also provide experimentally based correlations that 
are useful in defining parameters for the porous-media model; these correlations are 
documented in Gedeon Associates’ Sage Stirling-Code Manuals. These sources of 
experimentally based information were used to define the following terms and parameters 
needed in the non-equilibrium porous-media model: hydrodynamic dispersion, permeability, 
inertial coefficient, fluid effective thermal conductivity (including thermal dispersion and 
estimate of tortuosity effects), and fluid-solid heat transfer coefficient. Solid effective 
thermal conductivity (including the effect of tortuosity) was also estimated. Determination 
of the porous-media model parameters was based on planned use in a CFD model of 
Infinia’s Stirling Technology Demonstration Convertor (TDC), which uses a random-fiber 
regenerator matrix. The non-equilibrium porous-media model presented is considered to be 
an initial, or “draft,” model for possible incorporation in commercial CFD codes, with the 
expectation that the empirical parameters will likely need to be updated once resulting 
Stirling CFD model regenerator and engine results have been analyzed. The emphasis of 
the paper is on use of available data to define empirical parameters (and closure models) 
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needed in a thermal non-equilibrium porous-media model for Stirling regenerator 
simulation. Such a model has not yet been implemented by the authors or their associates. 
However, it is anticipated that a thermal non-equilibrium model such as that presented here, 
when incorporated in the CFD codes, will improve our ability to accurately model Stirling 
regenerators with CFD relative to current thermal-equilibrium porous-media models. 


Nomenclature 


A = heat transfer area (m 2 ) 

A S f = solid-fluid interfacial area (m 2 ) 
a = cross-sectional area (m 2 ) 

Cf = inertial coefficient 
c p = specific heat (J/kg-K) 
dh = hydraulic diameter (m) 
d w = wire diameter (m) 

/ = friction-factor 
h = heat transfer coefficient (J/m 2 -s-K) 

K = permeability (m 2 ) 
k = thermal conductivity (W/m-K) 
k tor = thermal tortuosity conduct, tensor (W/m-K) 
kdis = thermal dispersion conduct, tensor (W/m-K) 

L = length (m) 
m = mass flow rate (kg/s) 

Nu = h d h / k , Nusselt No. 
p = pressure (Pa) 

Pe = Re Pr, PecletNo. 

Pr = c p // / k , Prandtl No. 

Re = p u d h / u , Reynolds No. 
t = time (s) 

r h = hydraulic radius (m) = dh/4 
T = temperature (K) 

U = <u>, superficial, Darcian, or approach velocity 
u = average-flow-channel, or local, velocity 
V = volume (m 3 ) 

Greek Symbols 
P = porosity 
p = density (kg/ m 3 ) 
p = absolute viscosity (kg/m-s) 
v = dynamic viscosity (m 2 /s) 


I. Introduction 

The Advanced Stirling Convertor (ASC) is being developed as a light-weight, high-efficiency Advanced Stirling 
Radioisotope Generator (SRG) for use on potential NASA Space Science missions. Potential missions include 
providing spacecraft onboard-electric power for deep-space missions or power for lunar applications, or un-piloted 
Mars rovers. Advanced Stirling convertors would provide substantial performance and mass benefits for these 
missions and could also allow the use of Stirling radioisotope power systems for radioisotope electric propulsion. A 
combined Stirling convertor/cooler based on a high temperature heater head may enable an extended-duration Venus 
surface mission. Therefore, GRC is also developing advanced technology for Stirling convertors, aimed at 
improving specific power and efficiency of the convertor and the overall power system, and confirming and 
enhancing reliability. Performance and mass improvement goals have been established for these next generation 
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Stirling radioisotope power systems. Efforts are underway to achieve these goals, both in-house at GRC and via 
grants and contracts. These efforts include use of multi-dimensional (multi-D) Stirling computational-fluid- 
dynamics (CFD) models to enhance knowledge of component performance details and reliability, high-temperature 
materials, advanced controllers, low-vibration techniques, advanced regenerators, and a lightweight convertor. 12 

The objective of this paper is to report on definition of empirical parameters for an initial thermal non- 
equilibrium porous-media model intended for use in multi-D Stirling codes 3 " 5 for simulation of regenerators. 
Experimental data from a regenerator research grant 6 " 8 and experimentally based correlations derived from 
oscillating-flow rig test data 9 were used in defining parameters needed in the non-equilibrium porous-media model. 
The grant is led by Cleveland State University (CSU) with subcontractor assistance from the University of 
Minnesota (UMN), Gedeon Associates and Sunpower Inc. Determination of the particular porous media model 
parameters presented was based on planned use in a CFD model of Infinia’s Stirling Technology Demonstration 
Convertor (TDC). The thermal non-equilibrium porous-media model presented is considered to be an initial, or 
“draft,” model for possible incorporation in commercial CFD codes (which now contain thermal equilibrium porous- 
media models), with the expectation that empirical parameters or heuristic terms in the equations may need to be 
updated once resulting Stirling CFD model regenerator and engine results have been analyzed. The emphasis of the 
paper is on use of available data to define the empirical parameters (or closure models) needed in non-equilibrium 
porous-media models for Stirling regenerator simulation. Such a porous-media model has not yet been implemented 
in a Stirling CFD model by the authors or their associates. 

Two CFD codes, Fluent" and CFD-ACE* **, have been used at Glenn Research Center (GRC) for modeling of 
Infinia’s Stirling TDC. The porous-media models available in each of these codes are thermal equilibrium models, 
which assume that the solid matrix and the fluid are in local thermal equilibrium (thus at the same temperature) at 
each spatial interface between them within the porous medium. This is believed to be a poor assumption for the 
oscillating-flow environment within Stirling regenerators; some results supporting this belief are: (1) Experimental 
measurements in a large-scale replica of such a regenerator show solid and fluid temperatures differ by as much as 2 
degrees C over an engine cycle, as shown in Fig. 1 §§ . Characteristics of this large-scale replica and the tests were: 
90% porosity, stainless steel 304 welded screens, 200 layers of 6.3 mm x 6.3 mm mesh, 0.81 mm wire diameter; test 
section dimensionless parameters were similar to those of the TDC (the temperature measurements were made with 



Figure 1. Wire-Screen Regenerator (Considered to 
Approximate Random Fiber) and Test Results. 
The Temperature Measurements, shown at right, 
were made with a Hot Wire at the Axial Center of 
the 80 th Layer of Screen (from the Cold End). 
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hot-wire anemometry). (2) This same experimental data is shown in Figure 2, along with (a) a “microscopic” CFD 
simulation of the non-equilibrium temperatures of one layer of the test rig welded screen and adjacent fluid 10 , and 
(b) a CFD-ACE “macroscopic” equilibrium porous-media model simulation of the “lumped” solid plus fluid 
temperature over a cycle. Note that the variation of the calculated macroscopic thermal-equilibrium-model lumped 
solid/fluid temperature over the cycle (lying between the CFD microscopic and the measured solid temperatures in 
Fig. 2), is almost the same as that of the non-equilibrium solid temperature variations — due to the much greater 
specific heat of the solid. (3) Fig. 3 shows Sage 1-D TDC model predictions of the regenerator solid and fluid 


' 1 Internet website www.fluent.com or do internet search on “Fluent” 

** Internet website www.esi-group.com/SimulationSoftware/CFD-ACE/ or do internet search on “CFD-ACE” 

Simon, T., “Contributions to Regenerator Modeling,” Presentation for Regenerator Modeling Panel Discussion, 
2 nd International Energy Conversion Engineering Conference, Providence, RI, August 18, 2004. 


NASA/TM— 2006-2 14391 


3 




temperatures***. Differences in solid 
and fluid temperatures range up ~20 
degrees K at the ends of the 
regenerator and ~l-2 degrees at the 
center of the regenerator. These 
calculations, made by the first author at 
GRC, yielded an enthalpy flux loss 
over the cycle of about 22 W, or 
approximately 8% of the predicted 288 
W of heat into the hot-end of the 
convertor (predicted electrical power 
was about 56 W at ~83 Hz frequency); 
one would expect that neglect of the 
temperature difference between the 
regenerator solid and fluid, and the 
consequent entropy generation, would 
produce underestimation of 

Figure 2. CFD Microscopic Calculations of Solid and Fluid Temperatures, and regenerator enthalpy-flux losses. 

Macroscopic Equilibrium Solid+Fluid Temperature, deg. K— superimposed on Therefore these results, both 

the UMN Experimental Temperatures from Fig. 1 measurements and CFD simulations, 

suggest the assumption of solid- 

matrix/fluid thermal equilibrium would be inappropriate for the Stirling TDC regenerator. 


Ts-CFD Micro 
Tf CFD Micro 
Ts EXP 
Tf EXP 
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Regenerator Gas and Matrix Temperatures at 
Cold End 
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TDC Operating Conditions & Predictions 


• Operating Conditions (Not Design) were 

• Hot-End Temperature = 823 K (550 C) 

• Cold-End Temperature = 363 K (90 C) 

• Frequency of Operation = 82.5 Hz 

• Helium Mean Pressure = 2.59 MPa (376 psi) 

• Performance Predictions were 

• Engine Electrical Power = 55.9 We 

• Engine Efficiency = 19.4% 

• Heat Into Engine = 288 W 

• Regenerator Enthalpy Flux Loss = 22.4 W 

Regen. Enthalpy Flux Loss s j) 772 C ^ T gas dt 


Figure 3. Sage 1-D Thermal Non-Equilibrium Regenerator Model, Fluid & Solid Temperature Predictions for TDC 


A survey of the porous-media literature supports the need for thermal non-equilibrium porous-media models for 
thermal storage applications, as in Stirling regenerators. In Singh" it is noted that in a thermal storage system, heat 
transfer between the solid and fluid phases is a key phenomenon and both the phases are inherently in thermal non- 
equilibrium. Hence the corresponding models of heat transfer will necessarily be of the 2-energy-equation type, one 
equation each for the fluid and solid phase temperature. Results are shown in the paper of 1- and 2-energy-equation 
temperature calculations (i.e. using thermal- equilbrium and non-equilibrium models) for steady-flow through 


Tew, R., “Overview of Stirling Regenerator (Porous Media) CFD Modeling Problems,” Presentation for 
Workshop on Convective Heat Transfer in Porous Media, Hosted by NASA Glenn Research Center, North Olmsted, 
OH, April 28-29, 2004. 
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porous-media with thermal fields which evolve with time; 1- and 2-energy-equation temperature profiles are 
compared for glass/water and metal/air porous-media systems over ranges of Reynolds number and porous-media 
flow lengths. 

Amiri and Vafai 12 examined the validity of the thermal-equilibrium assumption via a one-dimensional study of 
incompressible flow through a packed bed. They noted that the transient response of the solid structure is crucial to 
the performance of packed beds when used for thermal storage applications, and the solid-to-fluid heat transfer 
coefficient is a major resistance to heat transfer at the interface between the fluid and the solid matrix. The 
“microscopic” governing equations were volume averaged over a representative elementary volume to arrive at the 
“macroscopic” transport equations, including separate fluid and solid phases; the resulting 1-D, incompressible, 
volume-averaged macroscopic equations used by Amiri and Vafai are a reduced form of the more general 
macroscopic equations given below in this paper. Temperature differences between fluid and solid were calculated 
for a range of Reynolds No.’s and Darcy No.’s (dimensionless permeability). It was found that the ratio of fluid-to- 
solid thermal conductivity was an important parameter in determining local temperature difference between fluid 
and solid, for a given set of operating conditions. This analysis was also for steady flow through the porous media 
(packed beds) with developing thermal fields. 

Minkowycz 13 analyzes departure from local thermal equilibrium in porous media due to a rapidly changing heat 
source. The analysis is done for a sinusoidally-varying imposed heat-flux at one of the boundaries for the cases of 
no-flow and steady-flow in the same direction as the imposed heat-flux. It is noted that local thermal equilibrium is 
an often-used hypothesis when studying heat transfer in porous media, but examination of non-equilibirum 
phenomena show that this hypothesis is usually not valid during rapid heating and cooling. Oscillatory plots of fluid 
and solid dimensionless temperatures as a function of dimensionless time over the cycle are shown for various 
combinations of dimensionless frequency, fluid and solid heat capacities, and dimensionless “Sparrow No.’s;” 
Sparrow No. is defined to be a product of one form of the more familiar Biot number and a dimensionless ratio of 
lengths. It is noted that in the presence of rapid surface heat input variations, local thermal equilibrium fails to exist 
and large errors can be realized, depending on the value of the thermophysical properties; and the value of the 
Sparrow number appears to be indicative of the presence, or not, of local thermal equilibrium for applications 
dealing with rapidly changing heat inputs. 

Muralidhar 14 analyzed the effect of pulsating flow of gas and the accompanying heat transfer within a 
regenerator made from wire-mesh screens, of the type used in Stirling cycle devices. A harmonic analysis was used 
to examine fully-developed but unsteady gas flow in a regenerator assuming a non-Darcy, thermally non- 
equilibrium porous medium. It is noted that one assumption often adopted in the published literature is the thermal 
equilibrium between fluid in the pores and the solid material of the porous medium; however, thermal performance 
of regenerators totally depends on heat transfer between fluid and solid phases. Therefore, the individual phases 
should have a local temperature difference from one to the other. Unfortunately, although temperature profiles and 
regenerator effectiveness are explored as functions of Reynolds No., pulsating frequency, thermal properties of gas 
and solid, and the mesh structure — the local temperature differences between gas and solid do not appear to be given 
in the paper. 

Gedeon 15 developed a 2-D CFD “Manifest” (for “manifold estimate”) computer code for simulation of Stirling 
regenerators and their manifolds (that was capable of extension to other Stirling components). In the process of 
developing this code, he derived a non-equilibrium macroscopic porous-media model, starting from the generalized 
Navier-Stokes equations. These integral-form porous-media equations were also converted to curvilinear form for 
use in Manifest. Gedeon used Manifest to study jetting into the porous matrix (as from adjacent heat exchangers) 
and the effect of side-inlet regenerator manifolds. Side-inlet manifolds appeared to greatly enhance the possibility 
that flow through the regenerator would be non-uniform and thus decrease regenerator performance. The form of 
Gedeon’s integral-form non-equlibrium porous-media equations differs from the differential form shown in this 
paper. 

Thermal equilibrium porous-media models, or 1 -energy-equation models, are consistent with the assumption of 
infinite heat-transfer coefficients between fluid and solid phases and completely neglect entropy generation due to 
the radial heat transfer between the fluid and solid phases (since the local temperature difference between the phases 
is zero). However, unless the solid heat capacity of the porous-medium is essentially infinite relative to that of the 
fluid, thermal-equilibrium-model hot-end and cold-end regenerator fluid exit temperatures can still vary over the 
cycle; so, in general, thermal equilibrium porous-media models of Stirling regenerators will still yield non-zero 
regenerator enthalpy flux losses over the cycle, though one would expect that these losses would be underestimated. 
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II. A Compressible-Flow Non-Equilibrium Porous-Media Model for CFD Codes 

Initial values of parameters needed for a macroscopic non-equilibrium porous-media model are defined below 
for use in modeling Technology Demonstration Convertor (TDC) regenerators in CFD codes. Fluid continuity, 
momentum and energy equations, and a solid-energy equation, are stated for reference in defining the parameters 
needed. Experimental and computational data generated under a regenerator research grant 6 ’ 8 , and data from an 
oscillating-flow test rig 9 are used in defining this initial set of parameters (i.e., in defining closure models). 

First define superficial (Darcian), approach, velocity as U = (») =ua / =w p . This is also the volume-flow 

rate through a unit cross-sectional area of the solid-plus-fluid. It is determined by averaging the velocity over a 
region that is small with respect to the macroscopic flow dimension but large with respect to the matrix-flow- 
channel size. In the above equation for superficial velocity, “ a ” is the average fluid-flow area within the total solid- 

plus-fluid cross-sectional area, a , , that the velocity is averaged over to determine the superficial velocity, U . The 
porosity is p . The average-flow-channel fluid velocity (or local porous-media velocity) is u , corresponding to the 
superficial velocity, t/=(w),and u = ( u)/p . 


A. Non-Equilibrium Porous-Media Conservation Equations 

The incompressible-flow, non-equilibrium, porous-media equations from Ibrahim and Rong 10 were used as a 
starting point; they used the volume-averaging approach described in Whitaker 16 in deriving their macroscopic 
porous-media equations. Rewritten for compressible flow 26 , and including a hydrodynamic dispersion term from 
Ayyaswamy" 1 ) the fluid continuity, momentum, and energy equations, respectively, written in terms of the 
superficial, or approach, or Darcian velocity are (where the brackets, < >, denote volume averaging): 


1 4<P> / < U >) i r , l 

b — + -pi v • k^> / <«x«>J 




(i) 


Ap) 

(mY 

K 


+ V 


P P 


(2) 


ttt 


(« )-{p) f ^k«>|<«) 
^ > - + jv * W («><*>'] = v • fc*v<r) / ]+ l)u ■ u 


(3) ss , 


§§§ **** 


V dV 


dA 


f 

V 


dV r 


where the wetted area per unit fluid volume, 
material derivative. Finally, the solid energy equation is: 


d(p\f 

(and r. = hydraulic radius) and — ,, — is a substantial or 


dt 


' ' ' Ayyaswamy, P. S., University of Pennsylvania, Private Communication, 2004. Ayyaswamy described this uu 
term as representing “the additional stresses due to the correlation between the velocity fluctuations associated with 
the hydrodynamic dispersion due to the solid matrix.” A closure model for the term is discussed later in this paper. 


‘ 1 ' The effective viscosity used in the “Brinkman term” of this momentum equation is different than the molecular 
viscosity, in general, as a result of the solid/fluid interfaces within the matrix. The Brinkman term computes 
momentum transport where there are velocity gradients within the flow. Examples would be the jetting of fluid from 
adjacent heat exchanger passages into the matrix, etc. 

sss d(p) p 

The pressure -work term, — , — , is the form for an ideal gas (helium, here). In general it should be written 


P T 


d{p) P 

dt 


For an ideal gas, the coefficient of thermal expansion, p 


1 

T ' 


Burmeister 17 suggests the form of the fluid energy viscous dissipation terms, 



(p) J 


£n 

Pk' 



here — which 


is consistent with the Darcy-Forchheimer terms in the momentum equation, (2). 
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( 4 ) 


d ( P s C s< r > S * ) 

dt 


r -| U/l /> 

V . fc.v(r) s J - h sf T) s -{T)f ) 


sf 1 1 

where the wetted area per unit solid volume, = — -t—* . 

r dV r, \—p 

s h r 

Kaviany 1 * * * has described a momentum equation very similar to equation (2) as a “semiheuristic” equation. That 
is, most of the terms are multi-dimensional terms based on derivation from fundamentals. However, the last two 
terms of equation (2) are based on a “heuristic” extension from a ID momentum equation (i.e., this is a “speculative 
extension” of an experimentally based ID equation to serve as a guide in study of multi-dimensional problems, that 
has not been experimentally verified for multi-dimensional problems). Thus, in this initial or “draft” non- 
equilibrium porous-media model, we will assume that permeability is isotropic (Even though we know that the wire- 
screen used in engines and the UMN experiments is not isotropic, and because of the method of fabrication of 
random fiber material used in modern engines, it is believed that most of the lengths of those fibers, also, lie in 
planes perpendicular to the flow direction). Thus, as this porous-media model is improved upon, it may be desirable 
in the future to use anisotropic permeability. There are forms of Equation (2) with terms similar to the last two in 
equation (2), but in which the permeability is a tensor quantity 16 . But, in equation (2), permeability is a scalar 
quantity assumed to be the same for all directions. Similar statements can be made for the inertial term multiplier of 


equation (2), i.e. C I 4k . 

The expressions or parameters in the above equations needing definition, for use in modeling Stirling 
regenerators, are: (1) the hydrodynamic dispersion term in the momentum equation; (2) scalar permeability, K , 

and (3) scalar inertial coefficient, Cy , in the momentum equation; (4) the fluid effective thermal conductivity, k p e , 


in the fluid-energy equation; (5) the solid effective thermal conductivity, k , in the solid-energy equation, and (6) 


the heat-transfer coefficient between fluid and solid, h ^ , appearing in both fluid and solid energy equations . Note 

that the thermal conductivity expressions in the energy equations are tensors, though they are expected to have non- 
zero values only on the diagonals. The products of velocity vectors in the momentum equation are also tensors. 


1 .Hydrodynamic Dispersion 

In the fluid momentum equation, equation (2) above, there is a hydrodynamic dispersion term: 

— (uu)suu 

p 


( 5 ) 


where ft is porosity; u is the average-channel fluid, or local, velocity inside the matrix, and u is the spatial-variation 
of the average-channel fluid velocity inside the matrix, uu is a tensor quantity which can be expanded in the form: 

uu = iiuu + iji7v + ikitw + jiv u + yjv v + jkv w + ki wu + kjvvv + kkvfff (6) 

For transport normal to the principal flow direction which is assumed to be the axial, i, direction, the mean 
velocity gradient is in the radial, j, direction and the term of interest (i.e. of significance) is u v , where v is the 
velocity component in that direction, = |u|j . Therefore the term from the above expression for transport normal to 
the flow is |«v| or v«| . The dispersion represented by this term is of axial momentum, pu , but this term is in 
the direction normal to u (in the direction of v’). 

In Niu 7 it was argued that this hydrodynamic dispersion could be equated to turbulent shear stress when the 
dispersion of interest is not in the direction of the mean flow, or («v) (Vv'^ (where the prime refers to temporal 


variations about the temporal average, the overbar refers to a temporal average, and < > refers to spatial average) 

— / — \ au 
and that (u'v'} could be modeled as £ m - ^h U where \uv/ - -s M and X = 0.02 ; anc j where d h = the 

hydraulic diameter, u = the in-matrix average velocity, and U=(u) is the superficial, Darcian, or approach velocity. 
Therefore, for use in equation (6), and in the momentum equation, (2): 


(8v) (vu) \ u ' v / 


P A P- 


-—0.02d. 17— = -—k d. U— 
h dr pi h dr 


pi M dr pi 


( 7 ) 
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2. Permeability and Inertial Coefficient 

In the fluid momentum equation, (2), the permeability, K , and the inertial coefficient, Cy , need evaluation for 

each type of porous medium. These coefficients can be evaluated for particular types of Stirling engine regenerator 
porous media via use of the friction-factor data from the oscillating-flow test rig given in Gedeon 9 under the 
assumption that the flow is quasi-steady 1 1 1 19 . The Darcy-Forchheimer steady-flow form of the 1-D fluid 
momentum equation, and a similar pressure-drop equation, but in terms of the Darcy friction factor, can be written, 
respectively: 

rWV (*) 

From the oscillating-flow test rig data, as summarized by correlations given in the Sage manuals 9 , it was determined 
that friction-factor correlations for random fiber and wire screen are of the form: 


Vp 

L 


fr 


a y 

+ 8 Re where Reynolds # , Re 


Pf ud h 


( 10 ) 


’D Re // 

And where, for random fiber, a=192, 5= 4.53, y= -0.067, and for woven screen, a=129, 8= 2.91, y= -0.103. 
Substituting the expression for friction-factor and the definition of Reynolds number from equations (10), into 
equation (9), and then equating the right-hand sides of equations (8) and (9), it can be determined that: 


K 


2 

a 


and Cy = 


8 Re' 




(ID 


A frequently used expression for hydraulic diameter of random fiber and wire screen, in terms of porosity and wire 
diameter is: 

P . 


1 -P 


( 12 ) 


For the welded stacked-screens used in the UMN regenerator test module 6 where /?= 0.9 and d = 0.81 mm or 


8.1E-4 m, then d, =7.29E-3 m. If it is also assumed that Reynolds No.’s are in the range from 25-100, as expected 

in the TDC regenerators, then equations (11) can be used to calculate values of permeability and inertial coefficient 
for TDC random fiber and the large screens used in the test module. The resulting TDC random fiber values are 
given in Table 1. The TDC random fiber results given in Table 1, in the rightmost column, are based on 
unidirectional flow tests of the entire TDC heater head (including heater, regenerator and cooler), taken from 
Wilson 19 The “UMN old” and “UMN new” values of permeability, in Table 1, are experimentally determined 
values of permeability and inertial coefficient, determined at an earlier and a later time®*. The “CSU Calcs.” are 
calculated values determined via microscopic CFD modeling of the UMN steady-flow test module*®. For the large 
scale-screens (see Table 1), results derived from the Sage empirical correlations and the microscopic CFD 
calculations are seen to be in relatively close agreement. The substantially greater differences of these values from 
the UMN test results have not been explained. The significance of differences such as these should be examined via 
a sensitivity study using a macroscopic porous-media model; this has not yet been done. 


3. Effective Fluid and Solid Thermal Conductivities 

In the fluid energy equation, (3), k = an effective fluid conductivity tensor, each element of which is, in 

J e 

general, a sum of components due to molecular conductivity, thermal tortuosity conductivity and thermal dispersion 
conductivity. It can be broken down into these components as follows: From equation (3), 

(13) 




\f. 


M 1 Tew, R. C., Private Communication, 2005. 

Simon, T., University of Minnesota, Private Communication, -2003. 

'® Ibrahim, M. B., Cleveland State University, Private Communication, 2003. 
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Table 1. Values of Permeability, K, and Inertial Coefficient, C f , for the Large-Scale UMN Wire Screen and TDC 
Random Fiber Regenerator Materials. 


Coefficient 

UMN Large-Scale Screens 
(d„,=8.1E-4 w) 

TDC Random Fiber 


UMN Old, 
Experimental 

UMN New, 
Experimental 

CSU Calcs. 

Sage Cor. 

Sage Cor. 

Unidirectional 
Flow Tests 

K(m 2 ) 

1.07E-7 

1.86E-7 

8.9E-7 

8.24E-7 

4.08E-10 

3.52E-10 

K/d w 2 

0.163 

0.283 

1.36 

1.26 

- 

- 

c f 

0.049 

0.052 

0.14 

0.13-0.11 

Re=25-100 

0.19-0.17 

Re=25-100 

0.154-0.095 

Re=25-100 


where: 

p ^Asf T dA = = k t or^( T Y 
f 

Equation (14) defines the thermal tortuosity thermal conductivity, k . 

j- 

-P f c p (f») , l dis . v(r)‘ 

Equation (15) defines the thermal dispersion thermal conductivity, . 

*fe * V ( T ) = F / * V ( r ) / + * tor * V ( r ) y + dis * V ( r ) y 


(14) 


(15) 


(16) 


Equation (16) defines the effective fluid thermal conductivity, k fe ■ 


and k r = k r + k r , + k = k r t + k 

je j J ,tor dis j ,stag dis 


(17) 


Where, in the above equation, the sum of the fluid molecular and thermal tortuosity conductivities are lumped 
together and called the fluid stagnant thermal conductivity. 

It is assumed that only the diagonal elements of the effective fluid conductivity tensor are non-zero. Then in 
terms of 3-D cylindrical coordinates which are appropriate for Stirling engine simulation, it is further assumed that: 



1 

> 

0 0 


k r +k r , +k J. 

j j ,tor,rr dis,rr 

0 0 

k r = 

fe 

0 

0 

k 0 

fe,99 

0 k f 

fe,xx J 


0 

0 

k f +k f,tor,99 +k dis,99 0 

0 k r -\-k r , ~^~k 

j j , tor, xx dis, xx _ 


k f ,stag ,rr + ^dis ,rr 

0 

0 


0 

^ f ,stag,3& + ^dis , 93 
0 


0 

0 

,stag , xx + k dis , xx 


(18) 


In the above tensor equation for fluid effective thermal conductivity, (18), molecular conductivity is isotropic, and 
fluid tortuosity and dispersion conductivities will be assumed anisotropic. In general, the tortuosity conductivity is 
different in different directions when the alternating fluid/solid geometry is different when looking in different 
directions. Therefore, when fluid molecular and thermal tortuosity conductivities are lumped together to form the 
stagnant thermal conductivity, the stagnant thermal conductivity is also, in general, different in different 
directions — as dictated by the geometry of the matrix. The thermal dispersion conductivity has contributions from 
an advective term and an eddy term. In the flow direction there would be both. Normal to the flow direction, there 
is only the eddy dispersion term. 

In the solid energy equation, (4), the effective solid conductivity, k ^ , is defined: 
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(19) 





/ \ 


k 0 0 

se,rr 


f\k,k 0 0 

J \s s,tor,rr ) 

/ \ 

k = 
se 

0 k se,99 0 

= 

0 ° 


0 0 k 

se,xx _ 


0 0 f{ k s\tor, XX )_ 


The molecular thermal conductivities can be assumed to be well known. That leaves to be determined: (1) fluid 
thermal dispersion and (2) fluid stagnant conductivities and solid effective conductivities — which are each functions 
of the molecular conductivity and the thermal tortuosity conductivity, for the fluid or solid; thus the form of the last 
matrix of equation (19). 

4. Thermal Dispersion Conductivity 

Niu 7 and other researchers, have measured eddy diffusivity (or transport) in, or very close to, porous media. If it 
is assumed that this is also equivalent to thermal dispersion due to eddies, as in Niu 7 , based on the Reynolds 
Analogy, then values for thermal dispersion due to eddies are given in Table 2, below, in terms of porous media 

hydraulic diameter, d h , and superficial, approach, or Darcian velocity, JJ =< u > . Although in general thermal 

dispersion by eddies is expected to be anisotropic, for this initial model it may be adequate to use the same 
relationship to calculate eddy dispersion in all directions. It should be noted, however, that in the flow direction, 
total dispersion is by eddies and by advection, where the advection term is dominant. For a direction normal to the 
flow, dispersion is only by eddies. To explore further possible differences of eddy dispersion in different directions 
refer to Niu " , McFadden and Simon . Once a macroscopic thermal non-equilibrium model is implemented in a 
Stirling CFD code, the significance of the differences in Table 2 can be explored via a sensitivity study. 


Table 2. A Comparison of Thermal Dispersion Coefficients from Several Methods. Taken from Niu 7 . 



Estimated Thermal Dispersion 

Porous media 

Current Direct 
Measurements at UMN, 
Niu 7 

e.. jj =*^=0.02 dJU or ^^=0.02 Pe 
M.eddy PfCp h kf 

Welded Screen 

Hunt and Tien 22 

k ^,yy =00011Pe 

kf 

Fibrous Media 

Metzger, Didierjean, 
and Maillet 23 

kjis.yy _( o.03 0.05 ).Pe and =0.073 Pe 1 59 

k f k f 

Packed Spheres 

Gedeon 9 

kdis.xj, 0.62 -2.91 ^dis.xx 

= 0.50 Pe P or a 0.06Pe for P = 0.9, Pe = 560 

k f 

Woven Screen 


5. Fluid-Stagnant and Solid-Effective Thermal Conductivity 

The fluid-stagnant and solid-effective thermal conductivities (each a function of the appropriate molecular and 
thermal tortuosity conductivities) are estimated based on the geometry of the matrix of interest. McFadden 20 
calculated a radial stagnant conductivity for the welded screens based on considerations of the geometry of the 
large-scale wire screens. Similar calculations are made for a fluid-saturated metal foam in Boomsma 24 . In the case 
of the Stirling TDC that has been modeled with CFD codes, a random fiber matrix is used in the regenerator. Since 
the details of the geometry for such a matrix are random, some assumptions must be made, as discussed below. 

For comparison with the calculations of McFadden 20 , mentioned above, the same molecular thermal 

conductivities used by McFadden will be used here (i.e. for air and stainless steel): So, assume that: k s =13.4 W/m- 
K for Stainless Steel 316, k j- = 0.026 W/m-K for air at standard temperature, and p =0.90. 

In the random fiber matrix most of the lengths of the fibers are believed to lie in planes perpendicular to the main 
flow axis. Therefore, initially make the assumption that for a 3-D CFD model, the effective solid+fluid conductivity 
in the radial and azimuthal directions (that would be appropriate for use in an equilibrium porous-media model) 
follows the parallel model defined below, and the effective solid+fluid conductivity in the axial direction follows the 
series model, also defined below. 

The parallel model for this lumped effective conductivity for fluid and solid, assuming all of the fibers run in the 
same direction (not including the fluid thermal dispersion) is: 
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k eff,s + f= k ff )+k S <i l -n (2°) 

Therefore, this effective solid+fluid conductivity for air and stainless steel combined, for a 90% porosity matrix, 
would be: 

-3 

k eff s+f = 126x10 w / mK ) (°-90) + (13.4 Wl mK) (0.1) 

= 0.0234 W / mK + 1.34 W / mK = 1.36 W / mK 


As already mentioned, the above would be an appropriate effective solid+fluid thermal conductivity for an 
equilibrium porous media model. However for a non-equilibrium porous media model, assume that the two terms 
on the right of the above equation represent the fluid stagnant thermal conductivity and the solid effective 
conductivity, respectively. That is, assume: 

-3 

k f stag = k f P = ( 26x10 w 1 mK ) (0.90) =0.0234 W / mK 


k se = k s (l— /0 = (1 3 .4 W / mK) (O.l) =1.34 W / mK 
However, in accordance with McFadden 20 : 


vv 


13.4/26x10 


-3 


515 and 


k eff, s+f lk f 


1.36/26x10 


52 


The value 52 is considerably larger than the 32.5 estimated by UMN for their wire screen in McFadden 20 , based on 
average geometrical considerations, rather than the parallel model. Since wire screen and random fiber are thought 
to have similar heat transfer properties, perhaps the solid part of the parallel model should be corrected to the UMN 
wire screen value by multiplication by the correction factor: 

kcoirect= 32.5/ 52=0.625 


Then the corrected parallel values for use for random fiber in the radial and azimuthal directions would be: 


k eff s+f = °- 022 4 + 1.36x0.625 = 0.873 (for equil., model) 
k =k (l-/l)o.625 = 1.34 W / mK x 0.625 = 0.838 (for non - equil. model, radial & azimuthal directions) 


For the axial direction, the lumped effective solid + fluid effective conductivity, k e jy s+ j , (which also doesn’t 

include thermal dispersion) should be substantially less than in the radial and azimuthal directions. Initially assume 
the series model mentioned in Ref. 14. That is: 


eff, s+f 


i 1 1 


p + { -p 

n 

^ 1 

V k j 

' f s' 

\26 


1 




OX 
-3 + 13.4 


0.0289 W / mK 


( 21 ) 


This axial effective solid+fluid thermal conductivity (not including thermal dispersion) is just slightly larger than the 
molecular fluid conductivity for air of 0.026 W/m-K and is probably too small since this series model assumes that 
the wires are not touching in the axial direction. 3-D CFD microscopic simulations of a representative elementary 
volume of the UMN regenerator by Rong , suggest that the value shown in equation (21) should be increased by 
a factor of 2.157. Therefore assume that: 


eff, s+f 



2.157 

1190 OX 
-3 + 1 3.4 


26x10 


(o.0289 IX / /« AT )(2 . 1 5 7 ) = 0.0623 If' / mK ( axial direction) (22) 


This resulting effective solid+fluid thermal conductivity in the axial direction would be appropriate for an 
equilibrium model. Seeing no obvious way to separate values for fluid and solid for a non-equilibrium macroscopic 
porous media model, based on the series model, we propose initially using this same value for both the fluid- 
stagnant and solid-effective conductivities in the axial direction in the non-equilibrium model, also — hoping that the 
overall axial effect might be reasonable. 

Recall that to get the total effective fluid conductivities in different directions, the thermal dispersion 
conductivity should be added to the fluid stagnant conductivities in the radial, azimuthal and axial directions. Also, 
for use in modeling the TDC regenerator, the above effective conductivities must be recalculated using the thermal 
conductivity of helium instead of air. 


Rong, W., Cleveland State University, Private Communication, 2005. 
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Literature searches for experimental effective fluid+solid effective thermal conductivities (appropriate for an 
equilibrium porous-media model) might prove useful. One attempt to locate such information via the website of the 
National Institute of Standards and Technology (NIST) was not successful. 


6. Heat Transfer Coefficient between Fluid and Solid Matrix Elements 

A good source of heat transfer correlations for wire screen and random fiber Stirling regenerator materials is 
Gedeon 9 . These correlations are based on experimental data from the NASA/Sunpower oscillating-flow test rig; this 
rig was designed and fabricated specifically for the purpose of determining friction-factor and heat-transfer 
correlations for use in Stirling device design and modeling. Heat transfer correlations in terms of Nusselt No., 
Peclet Number (Reynolds No. x Prandtl No.), and porosity are: 

For wire screen, Nu=(l. +0.99 Pe 066 ) /? 179 and for random fiber, A«=(L+1.16.Pe 0 - 66 )/? 2 ' 61 where 


Nu 


hd , 
l± 

k 


Pe = Re Pr = 


P ud h C P M 

p k 


Measurements by Niu 25 indicated that the above correlation and its application in a quasi-steady fashion is suitable 
for the portion of an oscillatory cycle when the acceleration is sufficiently weak or the flow is decelerating, but 
during strong acceleration, the unsteady measurements indicate a violation of the quasi-steady-flow assumption 
characterized by a lag between the heat flux from one phase of the matrix to the other and the temperature difference 
between phases. The measurements also showed an apparently poor mixing of the flow in the pore during the strong 
acceleration portions of the cycle. The Valensi number based on molecular (rather than eddy) viscosity) of the flow 
was 2.1, which is a bit higher than values in the general operating range of an engine regenerator (-0.23 for one 
engine). Thus, this unsteady effect may have been overestimated. The data also indicate that over the full cycle, an 
estimate given by the correlation above and applied as if the flow were quasi-steady is a reasonable approach for this 
initial or “draft” model. 


III. Concluding Remarks 

A set of transient, compressible-flow, conservation equations is summarized for reference in defining the 
parameters whose values are needed for a macroscopic, thermal-non-equilibrium porous-media model. Such a 
porous-media model is needed in existing commercial CFD codes (such as CFD-ACE and Fluent) in order to more 
accurately model the regenerator heat exchanger in Stirling engine devices (since only equilibrium porous-media 
models are now available in the Fluent and CFD-ACE codes). Available experimental information from large-scale 
wire-screen testing is used to define a hydrodynamic-dispersion term in the momentum equation. Experimental 
information is also used for definition of the permeability and inertial coefficients in the momentum equation and for 
the thermal-dispersion conductivity for the regenerator fluid. Methods are also outlined for estimating the stagnant- 
fluid and effective-solid thermal conductivities. Thus adequate information is presented for definition of an initial, 
or “draft,” non-equilibrium porous media model for use in CFD regenerator modeling of Stirling devices. It is 
anticipated that use of this initial model in CFD codes may demonstrate that further work on refinement of the non- 
equilibrium porous-media model and its parameters will be needed. 
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